********************************************************************************
********* Graphs for wedge heterogeneity by measures spending ******************


************************ for furnace
clear all
use "$maindir/Results/Model_Outputs/decompose_wedge.dta"

keep if substr(parm, -7, .)=="Furnace"
keep if substr(parm, 1, 1)!="0"

gen id_Furnace = _n

** merging with histogram data
merge 1:m id_Furnace using "$maindir/Results/Model_Outputs/histogram_data.dta"
drop if _n>12
drop _merge

** Holm correction for multiple hypothesis testing
sort p
multproc, puncor(0.05) pvalue(p) method(holm) critical(newp_holm) reject(reject_holm)
multproc, puncor(0.05) pvalue(p) method(bonferroni) critical(newp_bonfer) reject(reject_bonfer)
multproc, puncor(0.05) pvalue(p) method(hochberg) critical(newp_hoch) reject(reject_hoch)
multproc, puncor(0.05) pvalue(p) method(krieger) critical(newp_krieg) reject(reject_krieg)
multproc, puncor(0.05) pvalue(p) method(liu2) critical(newp_liu) reject(reject_liu)
multproc, puncor(0.05) pvalue(p) method(simes) critical(newp_simes) reject(reject_simes)


replace reject_simes = 0 if p==.


eclplot estimate min95 max95 id_Furnace, ///
    supby(reject_simes) ///
	estopts1(mcolor(gs4) msymbol(Th) msize(*2))  ///
	estopts2(mcolor(black) msymbol(S) msize(*2)) ///
	ciopts1(lcolor(gs4)) ///
	ciopts2(lcolor(black)) yline(0, lpattern(dash) lcolor(gray%50)) ///
	legend(off) ylabel(, glcolor(gray%10)) ///
	xtitle("Amount Spent (US$) on Furnace") ///
	graphregion(color(white)) bgcolor(white) ///
	ytitle("Performance Wedge (%)" "(compared to zero spending)") ///
	xlabel(1 "0" 2 "1-300" 3 "301-600" 4 "601-900" 5 "901-1200" 6 "1201-1500" ///
	7 "1501-1800" 8 "1801-2100" 9 "2101-2400" 10 "2401-2700" 11 "2701-3000" 12 ">3000" , ///
	labsize(small) angle(45)) yscale(range(-10 15) axis(1)) ///
	plot(bar freq_binFurnace id_Furnace, yaxis(2)	///
		fcolor(gray%25) fintensity(inten90) lwidth(none) ///
		barwidth(0.8) ytitle("Number of Homes", axis(2)) ///
		yscale(range(0 2000) axis(2)) ylabel(0(400)2000, axis(2)))	
graph export "$maindir/Results/Figures/wedge_furnace_overlay.pdf", replace
		
	
**** graph just to obtain the legend (used in all figures, later added in Latex)		
eclplot estimate min95 max95 id_Furnace, ///
    supby(reject_simes) ///
	estopts1(mcolor(gs4) msymbol(Th) msize(*2))  ///
	estopts2(mcolor(black) msymbol(S) msize(*2)) ///
	ciopts1(lcolor(gs4)) ///
	ciopts2(lcolor(black)) yline(0, lpattern(dash) lcolor(gray%50)) ///
	xtitle("Amount Spent (US$) on Furnace") ylabel(, glcolor(gray%10)) ///
	graphregion(color(white)) bgcolor(white) ///
	ytitle("Performance Wedge (%)" "(compared to zero spending)") ///
	xlabel(1 "0" 2 "1-300" 3 "301-600" 4 "601-900" 5 "901-1200" 6 "1201-1500" ///
	7 "1501-1800" 8 "1801-2100" 9 "2101-2400" 10 "2401-2700" 11 "2701-3000" 12 ">3000" , ///
	labsize(small) angle(45)) yscale(range(-10 15) axis(1)) ///
	legend(rows(3) order(2 "Non-Significant Coefficients" 4 "Significant Coefficients (after FDR corrections)" 5 "Number of Homes in Category")) ///
	plot(bar freq_binFurnace id_Furnace, yaxis(2)	///
		fcolor(gray%25) fintensity(inten90) lwidth(none) ///
		barwidth(0.8) ytitle("Number of Homes", axis(2)) ///
		yscale(range(0 2000) axis(2)) ylabel(0(400)2000, axis(2)))			
graph export "$maindir/Results/Figures/wedge_graph_legend.png", replace width(2500)
		


*********************** for air sealing
clear all
use "$maindir/Results/Model_Outputs/decompose_wedge.dta"

keep if substr(parm, -7, .)=="AirSeal"
keep if substr(parm, 1, 1)!="0"


gen id_AirSeal = _n

** merging with histogram data
merge 1:m id_AirSeal using "$maindir/Results/Model_Outputs/histogram_data.dta"
drop if _n>11
drop _merge

** Holm correction for multiple hypothesis testing
sort p
multproc, puncor(0.05) pvalue(p) method(holm) critical(newp_holm) reject(reject_holm)
multproc, puncor(0.05) pvalue(p) method(bonferroni) critical(newp_bonfer) reject(reject_bonfer)
multproc, puncor(0.05) pvalue(p) method(hochberg) critical(newp_hoch) reject(reject_hoch)
multproc, puncor(0.05) pvalue(p) method(krieger) critical(newp_krieg) reject(reject_krieg)
multproc, puncor(0.05) pvalue(p) method(liu2) critical(newp_liu) reject(reject_liu)
multproc, puncor(0.05) pvalue(p) method(simes) critical(newp_simes) reject(reject_simes)

replace reject_simes = 0 if p==.

eclplot estimate min95 max95 id_AirSeal, ///
    supby(reject_simes) ///
	estopts1(mcolor(gs4) msymbol(Th) msize(*2))  ///
	estopts2(mcolor(black) msymbol(S) msize(*2)) ///
	ciopts1(lcolor(gs4)) ///
	ciopts2(lcolor(black)) yline(0, lpattern(dash) lcolor(gray%50)) ///
	legend(off) ylabel(, glcolor(gray%10)) ///
	xtitle("Amount Spent (US$) on Air Sealing") ///
	graphregion(color(white)) bgcolor(white) ///
	ytitle("Performance Wedge (%)" "(compared to zero spending)") ///
	xlabel(1 "0" 2 "1-100" 3 "101-200" 4 "201-300" 5 "301-400" ///
	6 "401-500" 7 "501-600" 8 "601-700" 9 "701-800" 10 "801-900" 11 ">900", ///
	labsize(small) angle(45)) yscale(range(-7.5 5) axis(1)) ylabel(-7.5(2.5)5) ///
	plot(bar freq_binAirSeal id_AirSeal, yaxis(2)	///
		fcolor(gray%25) fintensity(inten90) lwidth(none) ///
		barwidth(0.8) ytitle("Number of Homes", axis(2)) ///
		yscale(range(0 1800) axis(2)) ylabel(0(360)1800, axis(2)))	
graph export "$maindir/Results/Figures/wedge_airseal_overlay.pdf", replace 



******************* for attic
clear all
use "$maindir/Results/Model_Outputs/decompose_wedge.dta"

keep if substr(parm, -5, .)=="Attic"
keep if substr(parm, 1, 1)!="0"

gen id_Attic = _n

** merging with histogram data
merge 1:m id_Attic using "$maindir/Results/Model_Outputs/histogram_data.dta"
drop if _n>10
drop _merge


** Holm correction for multiple hypothesis testing
sort p
multproc, puncor(0.05) pvalue(p) method(holm) critical(newp_holm) reject(reject_holm)
multproc, puncor(0.05) pvalue(p) method(bonferroni) critical(newp_bonfer) reject(reject_bonfer)
multproc, puncor(0.05) pvalue(p) method(hochberg) critical(newp_hoch) reject(reject_hoch)
multproc, puncor(0.05) pvalue(p) method(krieger) critical(newp_krieg) reject(reject_krieg)
multproc, puncor(0.05) pvalue(p) method(liu2) critical(newp_liu) reject(reject_liu)
multproc, puncor(0.05) pvalue(p) method(simes) critical(newp_simes) reject(reject_simes)

replace reject_simes = 0 if p==.
replace reject_krieg = 0 if p==.


eclplot estimate min95 max95 id_Attic, ///
    supby(reject_krieg) ///
	estopts1(mcolor(gs4) msymbol(Th) msize(*2))  ///
	estopts2(mcolor(black) msymbol(S) msize(*2)) ///
	ciopts1(lcolor(gs4)) ///
	ciopts2(lcolor(black)) yline(0, lpattern(dash) lcolor(gray%50)) ///
	legend(off) ylabel(, glcolor(gray%10)) ///
	xtitle("Amount Spent (US$) on Attic") ///
	graphregion(color(white)) bgcolor(white) ///
	ytitle("Performance Wedge (%)" "(compared to zero spending)")  ///
	xlabel(1 "0" 2 "1-300" 3 "301-600" 4 "601-900" 5 "901-1200" 6 "1201-1500" ///
	7 "1501-1800" 8 "1801-2100" 9 "2101-2400" 10 ">2400", ///
	labsize(small) angle(45)) yscale(range(-6 12) axis(1)) ylabel(-6(3)12, axis(1)) ///
	plot(bar freq_binAttic id_Attic, yaxis(2)	///
		fcolor(gray%25) fintensity(inten90) lwidth(none) ///
		barwidth(0.8) ytitle("Number of Homes", axis(2)) ///
		yscale(range(0 1200) axis(2)) ylabel(0(200)1200, axis(2)))			
graph export "$maindir/Results/Figures/wedge_attic_overlay.pdf", replace





******************* for wall insulation
clear all
use "$maindir/Results/Model_Outputs/decompose_wedge.dta"

keep if substr(parm, -7, .)=="WallIns"
keep if substr(parm, 1, 1)!="0"

gen id_WallIns = _n

** merging with histogram data
merge 1:m id_WallIns using "$maindir/Results/Model_Outputs/histogram_data.dta"
drop if _n>11
drop _merge
replace freq_binWallIns = 0 if id_WallIns==1

** Holm correction for multiple hypothesis testing
sort p
multproc, puncor(0.05) pvalue(p) method(holm) critical(newp_holm) reject(reject_holm)
multproc, puncor(0.05) pvalue(p) method(bonferroni) critical(newp_bonfer) reject(reject_bonfer)
multproc, puncor(0.05) pvalue(p) method(hochberg) critical(newp_hoch) reject(reject_hoch)
multproc, puncor(0.05) pvalue(p) method(krieger) critical(newp_krieg) reject(reject_krieg)
multproc, puncor(0.05) pvalue(p) method(liu2) critical(newp_liu) reject(reject_liu)
multproc, puncor(0.05) pvalue(p) method(simes) critical(newp_simes) reject(reject_simes)

replace reject_simes = 0 if p==.


eclplot estimate min95 max95 id_WallIns, ///
    supby(reject_simes) ///
	estopts1(mcolor(gs4) msymbol(Th) msize(*2))  ///
	estopts2(mcolor(black) msymbol(S) msize(*2)) ///
	ciopts1(lcolor(gs4)) ///
	ciopts2(lcolor(black)) yline(0, lpattern(dash) lcolor(gray%50)) ///
	legend(off) ylabel(, glcolor(gray%10)) ///
	xtitle("Amount Spent (US$) on Wall Insulation") ///
	graphregion(color(white)) bgcolor(white) ///
	ytitle("Performance Wedge (%)" "(compared to zero spending)") ///
	xlabel(1 "0" 2 "1-300" 3 "301-600" 4 "601-900" 5 "901-1200" 6 "1201-1500" ///
	7 "1501-1800" 8 "1801-2100" 9 "2101-2400" 10 "2401-2700" 11 ">2700" , ///
	labsize(small) angle(45)) yscale(range(0 24) axis(1)) ylabel(0(4)24, axis(1)) ///
	plot(bar freq_binWallIns id_WallIns, yaxis(2)	///
		fcolor(gray%25) fintensity(inten90) lwidth(none) ///
		barwidth(0.8) ytitle("Number of Homes" "(for 29.7% with non-zero spending)", axis(2)) ///
		yscale(range(0 420) axis(2)) ylabel(0(70)420, axis(2)))	
graph export "$maindir/Results/Figures/wedge_wallins_overlay.pdf", replace




********************* for windows
clear all
use "$maindir/Results/Model_Outputs/decompose_wedge.dta"

keep if substr(parm, -6, .)=="Window"
keep if substr(parm, 1, 1)!="0"

gen id_Window = _n

** merging with histogram data
merge 1:m id_Window using "$maindir/Results/Model_Outputs/histogram_data.dta"
drop if _n>12
drop _merge

** Holm correction for multiple hypothesis testing
sort p
multproc, puncor(0.05) pvalue(p) method(holm) critical(newp_holm) reject(reject_holm)
multproc, puncor(0.05) pvalue(p) method(bonferroni) critical(newp_bonfer) reject(reject_bonfer)
multproc, puncor(0.05) pvalue(p) method(hochberg) critical(newp_hoch) reject(reject_hoch)
multproc, puncor(0.05) pvalue(p) method(krieger) critical(newp_krieg) reject(reject_krieg)
multproc, puncor(0.05) pvalue(p) method(liu2) critical(newp_liu) reject(reject_liu)
multproc, puncor(0.05) pvalue(p) method(simes) critical(newp_simes) reject(reject_simes)

replace reject_simes = 0 if p==.

eclplot estimate min95 max95 id_Window, ///
    supby(reject_simes) ///
	estopts1(mcolor(gs4) msymbol(Th) msize(*2))  ///
	estopts2(mcolor(black) msymbol(S) msize(*2)) ///
	ciopts1(lcolor(gs4)) ///
	ciopts2(lcolor(black)) yline(0, lpattern(dash) lcolor(gray%50)) ///
	legend(off) ylabel(, glcolor(gray%10)) ///
	xtitle("Amount Spent (US$) on Windows") ///
	graphregion(color(white)) bgcolor(white) ///
	ytitle("Performance Wedge (%)" "(compared to zero spending)") ///
	xlabel(1 "0" 2 "1-200" 3 "201-400" 4 "401-600" 5 "601-800" 6 "801-1000" ///
	7 "1001-1200" 8 "1201-1400" 9 "1401-1600" 10 "1601-1800" 11 "1801-2000" 12 ">2000" , ///
	labsize(small) angle(45)) yscale(range(-3 12) axis(1)) ylabel(-3(3)12, axis(1)) ///
	plot(bar freq_binWindow id_Window, yaxis(2)	///
		fcolor(gray%25) fintensity(inten90) lwidth(none) ///
		barwidth(0.8) ytitle("Number of Homes", axis(2)) ///
		yscale(range(0 2000) axis(2)) ylabel(0(400)2000, axis(2)))	
graph export "$maindir/results/Figures/wedge_windows_overlay.pdf", replace 

